arXiv:1501.02808v2 [astro-ph.SR] 22Jun2015 


Draft version June 24, 2015 

Preprint typeset using style emulateapj v. 5/2/11 


SIGNATURES OF MRI-DRIVEN TURBULENCE IN PROTOPLANETARY DISKS: 

PREDICTIONS FOR ALMA OBSERVATIONS 

Jacob B. Simon^’^’^, A. Meredith Hughes'^, Kevin M. Flaherty'^, Xue-Ning BaP’®, and Philip J. Armitage^’’^ 

Draft version June 24, 2015 

ABSTRACT 

Spatially resolved observations of molecular line emission have the potential to yield unique con¬ 
straints on the nature of turbulence within protoplanetary disks. Using a combination of local non¬ 
ideal magnetohydrodynamic simulations and radiative transfer calculations, tailored to properties 
of the disk around HD 163296, we assess the ability of ALMA to detect turbulence driven by the 
magnetorotational instability (MRI). Our local simulations show that the MRI produces small-scale 
turbulent velocity fluctuations that increase in strength with height above the mid-plane. For a set 
of simulations at different disk radii, we fit a Maxell-Boltzmann distribution to the turbulent velocity 
and construct a turbulent broadening parameter as a function of radius and height. We input this 
broadening into radiative transfer calculations to quantify observational signatures of MRI-driven disk 
turbulence. We find that the ratio of the peak line flux to the flux at line center is a robust diagnostic 
of turbulence that is only mildly degenerate with systematic uncertainties in disk temperature. For 
the CO(3-2) line, which we expect to probe the most magnetically active slice of the disk column, 
variations in the predicted peak-to-trough ratio between our most and least turbulent models span a 
range of approximately 15%. Additional independent constraints can be derived from the morphology 
of spatially resolved line profiles, and we estimate the resolution required to detect turbulence on dif¬ 
ferent spatial scales. We discuss the role of lower optical depth molecular tracers, which trace regions 
closer to the disk mid-plane where velocities in MRI-driven models are systematically lower. 

Subject headings: accretion, accretion disks — (magnetohydrodynamics:) MHD — turbulence — 
protoplanetary disks 


1. INTRODUCTION 

Turbulence is the leading candidate for transport- 
ing angular momentum in gaseous accretion disks 
(jShaknra fc Sunvil^ 1197311. thus driving evolution o f 
the disk surface density (|Lvnden-Bell fc PringldlTOT^ l. 
In protoplanetary disks, turbulence is also important 
for planet formation processes due to coupling be¬ 
tween gas-phase turbulence and solid bodies. For 
small, aerodynamically coupled parti cles, turbulence can 
lead to radial and vertical di ffusion (jYoudin fc Lithwicld 
l2007t Iclarke fc PringldllQ^ . while simultaneously con- 
centrating particle s between vortices on small scales 
(jCnzzi et al.l 1200811. and within pressure maxima cre- 
ated by vortices (iBarre fc Sommerialll995ll or zonal flows 
(|Johansen et al.l[2009ll on large scales. For planets, whose 
coupling to the gas disk is gravitational rather than 
aerodynamic, turbulence can affe ct migration via the 
generation of turbulent t orques (jNelson fc Papaloizoul 
12004 iLubow fc Ida] 1201111 and i ts effect on the satu¬ 
ration of co-orbital reso nances (|Baruteau et al.l 120111 : 
iPaardekooper et ^1201111 . 
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Despite its fundamental importance, observational 
constraints on the nature of disk turbulence are limited 
and typically indi rect. In protoplanet ary disks, estimates 
of disk lifetimes (|Haisch et al.l [200lll a nd measurement 
of ag e -dependent stellar accr etion rates (|Gullbring et al.l 
1199^ iHartmann et al.l 1199811 can be interpreted as re¬ 
quiring a ~ 0.01 (adopting the a pr escription for turbu¬ 
lent angular momentum transport; iShaknra fc SnnvaevI 
Il973fl . The best prospects for improving upon this 
rough estimate lie in the detection of non-thermal turbu- 
lent broadening o f molecular lines, either i n the infrared 
(iCarr et al.l 1200411 or the submillimeter ([H ughes et al 


201lHCuilloteau et al.ll20l3:lde Greeorio-Monsalvo et al 


20131 1. These measurements are challenging. Very 
roughly, we expect velocity fluctuations to be of the or¬ 
der of 5v ~ O.lcs, where the sound speed Cg is 

itself smaller than the orbital velocity vk by a factor of 
the geometric thickness {h/r) ^ 0.05-0.15. Discerning 
the impact of turbulence on line profiles thus requires 
modeling of the kinematic and thermal structure of the 
disk at better than the 1% and 10% level, respectively. 
This is now feasible for a subset of well-resolved disk sys¬ 
tems with_accuratelyjMasured Keplerian rotation pro¬ 
files. iHughes et al.l (1201111 constrained turbulent velocity 
fluctuations in the upper layers of the disks surrounding 
TW Hydra and HD 163296 using the CO (3-2) line. For 
TW Hydra, they found an upper limit of 5v < O.lcg and 
in HD 163296 th e y mea sured 5v ^ O.dcg. Using CS lines, 
iCiiilloteau et al.l ( 201211 determined that 5v ^ 0.4 — 0.5cs 
in the molecular layer of the disk around DM Tau. Their 
measurements further suggested that this velocity is in¬ 
dependent of height above the mid-plane. 
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Observational determination of 5v yields direct insight 
into the magnitude of turbulence, but does not provide 
a model-independent constraint on angular momentum 
transport, particle concentration, or solid body migra¬ 
tion. The importance of turbulence in these processes 
depends upon the nature of the turbulence as well as 
its strength. Convection, for example, is a source of 
turbulenc e that is expected to be inefficient at trans - 
port fe.g.. lStone fc Balbuslll996t iLesur fc Ogilviell201Clll . 
while in the opposite direction, lar ge-scale internal mag¬ 
netic stresses or disk w i nds (e.g., Salmerone£_aljj2003 


Suzuki fc Inutsukall2009l: iBai fc Stonell2013aHLesur et al l 

2013t iFromang et al.ll2013ir can transport or extract an¬ 

gular momentum from a near-laminar disk. Even if we 
exclude winds, there are multiple possible drivers of gas- 
phase turbulence in protoplanetary disks, including grav¬ 
itational i nstability (iToomre 119641) . the vertical shear in¬ 
stability (INelson et al.l 12013 ), the sub critical baroclinic 
nstabilitv (|Lesur fc Papaloizoul 1201011 . and the magne - 


torotational instability IMRI: iBalbus fc HawlevI 1199811 . 
Discriminating between these alternatives may, however, 
be possible by comparing theoretical predictions for how 
6v varies with radius and height in the disk against ob¬ 
servational data. 

In this paper, we develop theoretical predictions for 
the radial and vertical variation of molecular line pro¬ 
files from a disk in which angular momentum transport 
is predominantly drive n by the MRI. Preliminary work 
bv ISimon et al.l (|2011af) has shown that the MRI results 


in turbulent velocity fluctuations Sv ^ O.lcs at high disk 
altitudes, and that the turbulent velocity is an increas¬ 
ing function of height above the mid-plane , particularly 
within the Ohmic dead zone (i Gammiell 19961). This is con - 
sistent with previous results (jFromang fc Nelsonll2006D . 
The trend in MRI-driven turbulence is distinctly different 
from that measured in simulations of gravitationally un- 
stable disks, in which Sv is almost indep endent of height 
(|Forgan et al.ll20ll: iShi fc Chiang|l20f^ . 

To make useful predictions for future observations re¬ 
quires substantial improvements over these prior calcu¬ 
lations. First, we need a more accurate treatment of 
non-ideal magnetohydrodynamic (MHD) effects, which 
are central t o the outcome of the MRI in protoplan¬ 
etary disks (lArmitagd 120111) due to the l ow ioniza - 
tion fraction in these systems (iWardlel [2007t iBail 120111) . 
At high densities. Ohmic diffusion can prevent cou¬ 
pling of the magnetic field t o the dis k gas, leading to 
the dead zone paradigm of IGammid (119961) in which 
MRI-turbulent surface layers are ionized due to X- 
rays and cosmic rays and surround a low-ionization 
dead zone in which the re is only very weak accretion 
(|Fleming fc Stonel 120031) . At low gas densities, ambipo- 
lar diffusion reduces the efficiency of the MRI due to 
weak collisional coupling between ions and neutrals , lead¬ 
ing to the ambipolar damping zone (IBai fc Stoiid 120111 : 
iMohantv et al.l[2M^ iSimon et al.ll2013bllaD in the outer 
disk, and qu enched turbulence in the low density d isk 
atmosphere (jBai fc Stonel I2013bl : ISimon et al.ll20i3bl R). 
At intermediate densities, the Hall effect leads to a qual¬ 
itatively different behavior in which the magnitude and 
nature of magnetohydrodynamic transport depends on 
the orientation of the net vertical field with respect to 
the disk angular momentum v ector (|Kunz fc Lesui][2013t 
iLesur et al.l 12014 lBaill2014a|) . Here, we calculate a se¬ 


ries of local disk models, incorporating both Ohmic and 
ambipolar diffusion, whose physical properties (density, 
temperature, accretion rate) are tailored to match spe¬ 
cific radii in the disk around HD 163296 (at the radii 
of greatest interest observation ally, we bel ieve that the 
Hall effect is sub-dominant; e.g., [Bai2014bL Simon 2015, 
in prep). Second, we need to derive the directly observ¬ 
able quantities (molecular line profiles and channel maps) 
from the primary theoretical output (velocity fields at 
different levels within the disk). To carry out this task, 
we use the LIME radiative transfer code in conjunction 
with an empirical model for the vertical thermal struc¬ 
ture of the disk to generate molecular lines and channel 
maps. This approach allows us to make the first real¬ 
istic predictions for sub-mm observations of HD 163296 
under the assumption that any turbulence in the system 
derives from the MRI. These predictions will be testable 
with forthcoming high resolution, high sensitivity ALMA 
observations of the HD 163296 disk. 

The outline of this paper is as follows. In Section [2l 
we describe the empirical model for the HD 163296 disk 
that we use throughout this work. In Section [31 we first 
explain our numerical methodology and then discuss the 
results from our simulations, including the extraction of 
the turbulent velocity for input into the LIME code. Our 
radiative transfer calculations and observational predic¬ 
tions are presented in Section |4l and we summarize and 
present our conclusions in Section |S] 

2. DISK MODEL 

Our numerical simulations and radiative transfer cal¬ 
culations are based on a representative model for the disk 
around HD 163296, which is broa dly consistent with pre - 
vious observations of this source ([Rosenfeld et al.ll2013D . 
The disk model is needed for two distinct purposes in 
this work. First, the surface density and mid-plane tem¬ 
perature (which sets the vertical scale height) at a small 
number of discrete radii serve as input to the local MHD 
simulations that we use to determine the predicted turbu¬ 
lent velocity. The simulations are isothermal, and hence 
are characterized by a single temperature and have an 
approximately Gaussian vertical density profile (modi¬ 
fied only by the effects of magnetic pressure). Second, 
an axisymmetric disk model is needed as input to the 
radiative transfer calculations that we use to determine 
the predicted line shapes and channel maps. A verti¬ 
cally isothermal disk structure would yield results that 
are grossly inconsistent with observations because we are 
interested in line emission that originates from near the 
surface of the disk where the temperature is substantially 
elevated. Accordingly, for the radiative transfer calcula¬ 
tions we use a version of the disk model in which T is a 
function of height as well as radius, though we continue 
to assume a Gaussian density profile. This approach is 
computationally expedient - non-isothermal disk simu¬ 
lations are significantly more difficult to develop and run 
- but it carries a cost in terms of consistency; neither 
the radiative transfer model nor the matching between 
the simulations and the radiative transfer can be fully 
self-consistent. 

The vertical temperature structure is consistent with 
warm surface layers surrounding a cooler mid-plane. We 
first define the “mid-plane” and “surface layer” temper¬ 
atures, respectively, as 
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"^”’°(l55Au) 

Ts{r) = Ts,o (200AU) 
The complete gas temperature is 


r T(r,z} = Ts(r) 

where Zq and 5 depend on the radius r via, 


-(s5au)' 


and 


6 = 0.0034 


(®) 


- 200 


+ 2.5. 


( 1 ) 


TABLE 1 

Model Parameters 


( 2 ) 


ld>. 

\z\ < . 

(3) 


Parameter 

Value 

7m ,0 

19 K 

Ts,0 

40 K 

+tdisk 

0.09 Mq 

-^star 

2.3 Mq 

ID 

0.8 

Zq Rc 

115 AU 

Zq 

2.37 


of z because the gas density is a simple Gaussian in z. 
Thus, after some basic arithmetic, we set the CO abun¬ 
dance to this low value if Erfc[z/77(r)] < 10^^mH/E(r) 
(4) and if Erfc[—z/i?(r)] < 10^^mH/E(r) (here, the pre fac¬ 
tor is assumed to have appropriate units to make the 
ratio dimensionless). 

Our model disk parameters are presented in Table [TJ 


( 5 ) 


3. LOCAL DISK SIMULATIONS 


The surface density for HD 163296 is well fit by a power 
law with an exponential drop off at large radii: 


E(r) 


Afdisk 


2 — 7_d 

27ri?2 




( 6 ) 


from which we can calculate the volume mass density. 


p(r,z) 


S(r) 

^/nH{r) 


exp 



where the vertical scale height H{r) is 


( 7 ) 


H = 


'/2cs 

fl 


( 8 ) 


The sound speed (which we take to be isothermal here 
given the Gaussian density profile) is 


Cs 


^boltz7ni(r ) 




( 9 ) 


where Em is the mid-plane temperature profile from 
Equation o, and 


appropriate for Keplerian rotation, fcboitz is Boltzmann’s 
constant, ttih is the mass of hydrogen, and G is the grav¬ 
itational constant. 

Einally, our radiative transfer model includes the abun¬ 
dance of CO relative to hydrogen, Xqoj in order to pro¬ 
duce synthetic observables using CO lines. In most re¬ 
gions, we assume Xco = 10“^. If the temperature is less 
than 19K, we set Xqo = 10“^^, representing freeze-out 
of CO onto dust grains. Photo-dissociation also depletes 
CO in the surface layers. When N < 5 x 10^° cm“^ 
(where N is the integrated number density from the sur¬ 
faces of the disk), we set Xco = 10“^^. We can translate 
this to a condition on the complementary error function 


Our numerical simulations are a series of local, disk 
simulations of an accretion disk patch placed at sev¬ 
eral different radii in the model for HD 163296 described 
above. Here, we describe the numerical algorithms and 
initial conditions for these simulations. In this paper, 
there will be two sets of units. The first set of units will 
be code units, which we will refer to when we are talking 
about specifics of numerical simulations. The second set 
of units will be physical units, which will apply to our ra¬ 
diative transfer calculations and all other analysis from 
that point on. 


3.1. Numerical Method and Setup 

Our simulations use Athena, a second-order accurate 
Godunov flux-conservative code for solving the equations 
of MHD. Athena uses the dimensionally unsplit corner 
transport upwind method of lColell^ (II990I1 coupled with 
the third-order in space pi ecewise parabolic method of 
iColella fc WoodwardI |I984 ) and a constrained transport 
('CT: lEvans fc Hawlevlll9^ ) algorithm for preserving the 
V ■ B — 0 constraint. We use the HE LD Riemann solver 
to calculate the nu merical fluxes (jMivoshi fc Kusanol 
120051: lMignonell200^ . A detailed description of the base 
Athena algo rithm and the results of v a rious test problems 
Siie giv en in Gardiner fc Stone (|2005ll . [Gardiner &: Stonel 
(|200k) . and Stone et al.l ([2008 1. 

The simulations employ a local shearing box approx¬ 
imation. The shearing box models a co-rotating disk 
patch whose size is small compared to the radial distance 
from the central object, Rq. This allows the construc¬ 
tion of a local Cartesian frame (x, y, z) that is defined in 
terms of the disk’s cylindrical co-ordinates {R, (j), z') via 
X = {R — Rq), y = Rocj), and z = z'. The local patch 
co-rotates with an angular velocity D corresponding to 
the orbital fr e quenc y at Rq, the center of the box; see 
iHawlev et al.l (1199511 . The equations to solve are: 


dp 


+ V • (pv) = 0, 


^ + V • {pvv - BB) + v(p+ 


= 2qpirx — pfrz — 2S1 X pv 


(11) 

( 12 ) 
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dB 


V X (t> X B) = -V X (?7 aJj_ + voJ), 


(13) 


where p is the mass density, pv is the momentum density, 
B is the magnetic field, P is the gas pressure, and q is 
the shear parameter, defined asq = —dlnQ/dlnR. We use 
q = 3/2, appropriate for a Keplerian disk. For simplic¬ 
ity and numerical convenience, we assume an isothermal 
equation of state P = pc^, where Cg is the isothermal 
sound speed. From left to right, the source terms in 
equation m correspond to radial tidal forces (gravity 
and centrifugal), vertical gravity, and the Coriolis force. 
The first electromotive force (EMF) term on the right- 
hand-side of equation m is ambipolar diffusion, which 
consists of diffusivity tja = multiplied by the 

component of the current density J perpendicular to the 
magnetic field, Jj_. Here, pi is the ion density, and 7 is 
the coefficient of momentum transfer for ion-neutral col¬ 
lisions. Ohmic diffusion is included via the second EMF 
term on the right-hand-side and is the diffusivity po mul¬ 
tiplied by J. Our system of code units has the magnetic 
permeability equal to unity, and the current density is 


J = V X B. 


(14) 


Numerical algorithms f or integrating thes e equa tions 
are described in det ail in IStone fc Gardineil (1201011 (see 
also the Appendix of ISimon et al.ll2011bll . The y bound¬ 
ary conditions are strictly per iodic, whereas thex bound¬ 
aries are shearing periodic (iHawlev et al.l The 

vertical boundary con ditions are the mod ified outflow 
boundaries described in lSimon et all (|2013a|l . The EMFs 
at the radial boundaries are properly remapped to guar¬ 
antee that the net vertical m agnetic flux is conserved to 
machine p recision using CT (IStone fc Gardineill20f^. 


As in iBai &: Stone! (1201111. Simon et al. (1201 lall 


iSimon et al.l (|2013bll . and iSimon et al.l (|2013all , both 
Ohmic and ambipolar diffusion are implemented in a 
first-order in time operator-split manner using CT to pre¬ 
serve the divergence free condition with an additional 
step of remapping Jy at radial shearing-box bound- 
aries. The super ti me-stepping (STS) technique of 
lAlexiades et al.l (|1996ll has been implemen ted to accel¬ 
erate o ur calculations (see the Appendix of iSimon et al.l 

(HoTsB)). 

We calculate the Ohmic and ambipolar diffusivities 
by interpolating from pre-computed diffusivity tables 
based on equilibri um chemistry , follow ing the method¬ 
ology described in IBai fc Stone! (I2013b(l and subsequent 
works (|Baill2013L l2014aH . The diffusivity tables give po 
and pa/B^ as a function of density and ionization rate 
at fixed temperature. They are independent of B for 
the regimes we consider in this paper (in the absence 
of abundant small gr ains). The chemi c al rea ction net¬ 
work is described in IBai fc GoodmanI (1200911 and iBail 
dsnni), and recently undated in IBail (l2Q14all with the 
latest version of the UMIST database (|McElrov et al.l 
1201311 . We consider both grain-free chemistry, as well as 
a chemistry model containing 0 . 1 /im grains with abun¬ 
dance of 10 “'*. In the simulations, the ionization rates 
are obtained based on the horizontally averaged col¬ 
umn density to the disk surface, E. We include con- 
tribut ions from stellar X-ray s using standard prescrip¬ 
tions (llgea fc Glassgoldlll999ll. w ith fitting formulas pro¬ 
vided by^a^^^oodma^ 200 ^ , and cosmic rays with 


ionization rate ^cr = 10“*”^exp [—E/(96 g cm“^)]s“* 
(jUmebavashi fc Nakanol 1198111 . For X-ray ionization, 
we use an X-ray luminosity Lx = 4 x 10^® erg s“* 
and X-ray temperature correspondi ng to 1 keV, s pecif¬ 
ically appropriate to H D163296 (|Swartz et al.l 120051 : 
iGiinther fc Schmittl 1200911 . In the surface layer, we fur¬ 
ther include the effect of far UV ionizat i on ba sed on 
calculations from fPerez-Becker fc Chiangl (120111 1. The 
FUV ionization is assumed to have a constant penetra¬ 
tion depth of Epuv in the range of 0.01 — 0.1 g cm“^. 

The degree of FUV ionization is critically important 
for MHD turbulence in the outer disk, because only FUV 
photons can yield ionization levels in the upper disk lay¬ 
ers that are high enough to prevent MRl quenching by 
ambipolar diffusion. We treat the effect of FUV ioniza- 
ti on on ambipo l ar diffu sion using the same procedure as 
in ISimon et al.l (|2013all . We define the ambipolar diffu¬ 
sion Elsasser number. 


Am.f, (15) 

which corresponds to the number of times a neutral 
molecule collides with the ions in a dynamical time 
(r2“*). Since in the FUV layer. Am ^ 1 (i.e., ambipolar 
diffusion is weak) at a l l radii, we simply adopt the func¬ 
tion from ISimon et al] (l2013all . which is sufficient for our 
purposes. In the FUV ionization layer. Am is 


Ampuv 



The Am value below the FUV ionization layer, Ammid is 
determined by the diffusivity table, as described above. 
Piecing these two regions together, we have the following 
complete functional form for Am, 


Am = 


Ampuv 

Auirnid + |AmFuvS'+(2:) 
< AlUmid 

AlUmid + |AmFUvS'“(2 ) 
_ Aiufuv 


z > zt + Az 

zt — nAz < z < Zt + Az 
Zb + nAz < z < Zt — nAz 
Zb — Az < z < Zb + nAz 
z < Zb — Az 


(17) 

where Zt and Zb are the top and bottom layers of the FUV 
ionization layer, respectively, and are calculated by inte¬ 
grating at each time step the horizontally averaged col¬ 
umn density from the boundary towards the mid-plane 
until Epuv is reached. 5 ''''( 2 ) and S~{z) are smoothing 
functions defined as 


g+(^)^l + Erf aT* )’ 

5“ {z) = l- Erf , (19) 

Here, n = 8 and Az ranges from O.IH to 0.05H, de¬ 
pending on the simulation. These numbers were chosen 
to give a reasonably resolved transition region between 
Am = Ammid and Ampuv- 

We center the shearing box simulations on several radii 
in the mid-to-outer regions of HD 16329 6, at 10 AU, 
30 AU, and 100 AU. Rec ent results (e.g., IBai fc Stonel 
l2013al: ISimon et al.ll2013bl ial) point to the importance of 
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TABLE 2 

Shearing Box Simulations 


Label 

Radius 

do 

Spuv 

Grains? 

Q 

M 

Simulation Category 


(AU) 


(g cm-2) 



{Mq/ yr) 


R10-B4-FUV0.1-NG 

10 

10^ 

0.1 

No 

0.0071 

2 x 10“'^ 

Non-dynamo 

R10-B5-FUV0.01-NG 

10 

10® 

0.01 

No 

0.00081 

2 x 10“® 

Non-dynamo 

R10-B5-FUV0.1-NG 

10 

10® 

0.1 

No 

0.0024 

2 x 10“® 

Dynamo 

R10-B5-FUV0.01-G 

10 

10® 

0.01 

Yes 

0.00078 

2 x 10“® 

Non-dynamo 

R30-B4-FUV0.01-NG 

30 

10^ 

0.01 

No 

0.0052 

1 x 10“'^ 

Non-dynamo 

R30-B4-FUV0.1-NG 

30 

10^ 

0.1 

No 

0.0099 

2 x 10“'^ 

Non-dynamo 

R30-B5-FUV0.01-NG 

30 

10® 

0.01 

No 

0.0013 

2 X 10“® 

Non-dynamo 

R30-B5-FUV0.1-NG 

30 

10® 

0.1 

No 

0.0032 

2 X 10“® 

Dynamo 

R30-B5-FUV0.01-G 

30 

10® 

0.01 

Yes 

0.0013 

2 X 10“® 

Non-dynamo 

R100-B4-FUV0.01-NG 

100 

10^ 

0.01 

No 

0.0075 

9 X 10“® 

Non-dynamo 

R100-B4-FUV0.1-NG 

100 

10^ 

0.1 

No 

0.027 

1 X 10“’^ 

Dynamo 

R100-B5-FUV0.01-NG 

100 

10® 

0.01 

No 

0.0026 

9 X lO”'^ 

Dynamo 

R100-B5-FUV0.1-NG 

100 

10® 

0.1 

No 

0.0044 

1 X 10“® 

Dynamo 

R100-B5-FUV0.01-G 

100 

10® 

0.01 

Yes 

0.0025 

8 X 10“® 

Dynamo 


the strength of the net vertical magnetic field thread¬ 
ing the domain as well as the depth to which FUV pho¬ 
tons ionize the upper disk layers. As such, we choose 
these parameters to explore in our parameter study. In 
particular, we examine initial field strengths, defined by 
the gas to magnetic pressure ratio at the mid-plane, 
/3o, of /3o = 10'^ and 10®. We explore the limiting 
cases for the column to which FUV photons penetrate , 
taken from the results of IPerez-Becker fc Chian^ (|2011f) : 
Spuv = 0.01 g cm“^ and Spuv = 0.1 g cm“^. We also 
run one case at each radius that includes the effects of 
grain chemistry on the diffusivity coefficients, though we 
don’t find a significant difference between this simulation 
and the corresponding simulations with no grains. 

Aside from the initial field strength and details of the 
diffusion profile that depend on the inclusion of grain 
chemistry and the depth of FUV ionization, all simu¬ 
lations start from the same initial conditions. The gas 
density is in hydrostatic equilibrium for an isothermal 
gas, 

p{x, y, z) = poexp ^ ( 20 ) 

where po = 1 is the mid-plane density in code units, and 
H is the scale height in the disk defined by the mid-plane 
gas temperature (see Equations m and ®) 

In code units, the isothermal sound speed, Cs = 7.07 x 
10 “^, corresponding to an initial value for the mid-plane 
gas pressure of Po = 5 x 10“^. With — 0.001, the value 
for the scale height is i7 = 1. A density floor is imposed 
as too small a density leads to a large Alfven speed and 
a very small time step. The value for the density floor 
depends on the simulation but ranges from 10 “® to 10 “'* 
of the initial mid-plane value. 

In all runs, the initial magnetic field is a net vertical 
field. It is well known that for purely vertical field, the 
MRI sets in from a transient channel flow. For relatively 

I 


strong net vertical magnetic flux, the channel flow is so 
strong as to cause numeri cal problems and/or disk dis¬ 
ruption in the simulations (jMiller fc Stonel[200r)fl . To cir¬ 
cumvent such potential difficulties, we add a sinusoidally 
varying vertical field component so that, 


r. 1 • / 

^27r 

1 + -sin 

—x 

2 

) 


where is the domain size in the x dimension, and 

( 22 ) 

Here, Bq is the net vertical magnetic field, characterized 
by /3o, described above. With the asymmetry introduced 
by the extra sinusoidal variation in t he vertical field, the 
strong growth of cha nnel flows (see iHawlev et al.l 119951 : 
iMiller fc Stond [^OOClO is suppressed at early stages, and 
the simulation can integrate beyond the initial transient 
without numerical problems. All other magnetic field 
components are initialized to be zero. 

To seed the MRI, random perturbations are added to 
the density and velocity components. The amplitude of 
these perturbations are 6p/po= 0.01 and Svi = O.OOdcs 
for i = x,y^ z. All simulations are carried out at a resolu¬ 
tion of 36 zones per H and at a box size of AH x SH x SH. 

To summarize, we have three free parameters: the ini¬ 
tial magnetic field strength, the column to which the 
FUV photons penetrate, and the inclusion of dust grains 
in our chemistry calculation of the diffusivity profiles. 
The runs are listed in Table [H The label for each run 
is given by its radial location, field strength, the FUV 
depth, and the chemistry calculation. So, for example, 
run R10-B4-FUV0.01-NG is placed at a radius of 10 AU, 
has a magnetic field strength defined by /Iq = 10 "*, has 
an FUV penetration column of 0.01 g cm“^, and has no 
grains “NG”. 
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3.2. Simulation Results 

In this section, we first discuss the basic properties 
of our simulations, and how they compare to existing 
results in the literature. We then describe our method¬ 
ology for extracting a mean vertical profile of turbulent 
velocity from our runs for subsequent use in the radiative 
transfer calculation of molecular line profiles and other 
observables. 


3.2.1. Preliminaries 

To set the stage for our turbulent velocity character¬ 
ization below, we first examine our simulations through 
sever al standard diagnostics. Th e first is the a parame¬ 
ter of iShakura fc SunvaevI (119731 1. calculated as the time 
average of the density weighted R(j) {xy in shearing box 
coordinates) component of the stress tensor. 


a = 


{pVxSv, 


y BxBy) 


{pci) 


(23) 


where the angled brackets denote a volume average and 
the over-bar denotes a time average. The value of a 
for each simulation is given in Table [2] and generally lies 
within the range ^ 0.001 — 0.01, consist ent with previous 
shearing box stud ies of the MRI ()e.g.. lSano et al.ir2004 
iSimon et al.l 1201211. The a values also agree with the 
simulations of ISimon et al.l (l2013all . which explored the 
nature of MRI turbulence in the presence of ambipolar 
diffusion and a net vertical field in the outer regions of a 
disk with a Minimum Mass Solar Nebula density profile. 

The simulations can be classified into two categories 
based on the evolution and structure of the magnetic 
held. The hrst type of solution, which we refer to as 
a “dynamo” solution, shows the standard MRI dynamo 
behavior for the magnetic held. In particular, the mean 
radial and toroidal helds hip sign with a period of ^10 or¬ 
bits in the active regions that surround the Ohmic and/or 
ambipolar diffusion dominated mid-plane region. The 
second type of solution, which we call “non-dynamo”, 
consists of a mean magnetic held that is approximately 
constant in time. 

Figure [T] shows the space-time evolution of the horizon¬ 
tally averaged toroidal held for two non-dynamo simula¬ 
tions (top two panels) and one dynamo simulation ( bot- 
tom panel). In previous work ( Simon et al.l[^13all we 
described the non-dynamo simulations as being “quasi- 
laminar”, but here we choose a different naming conven¬ 
tion to emphasize that these solutions are not necessarily 
laminar in nature. In particular, for most of our non¬ 
dynamo simulations, there is a signihcant Maxwell stress 
contribution from scales other than the largest scale in 
the box. This is shown in Figure [2l which shows the 
time and horizontally averaged vertical prohle for the xy 
component of the Maxwell stress, {—BxBy), for three 
separate simulations. Also plotted in the figure is the 
small scale contribution to the Maxwell stress, calculated 
by subtracting the Maxwell stress at the largest scale, 
— {Bx){By), from {—BxBy). The degree to which these 
stresses overlap represents how much of the stress resides 
at large scales as opposed to turbulent fluctuations. The 
left and middle panel show the two non-dynamo simula¬ 
tions from Fig. [1] whereas the panel on the right shows 
the dynamo simulation. In the non-dynamo cases, the 


Maxwell stress can be largely dominated by large scale 
stress (as in the left panel), but in general these simula¬ 
tions have some small scale turbulence towards the mid¬ 
plane region (middle panel). In the dynamo simulation, 
the Maxwell stress is generally large scale at large \z\, but 
has a much greater contribution from small scale stress 
compared to the non-dynamo simulations. Furthermore, 
the mid-plane regions of the dynamo simulations are al¬ 
ways dominated by small scale turbulent stress. The 
simulation category is displayed in Table [2J 

One feature of note is the strong asymmetry in the 
small-scale turbulent stress of the left panel of Fig. [5J 
We are not entirely sure why this asymmetry exists, but 
we believe that it results from a combination of stochas¬ 
tic variation in the magnetic field strength, which gen¬ 
erally plays a role in driving small scale stress, and rela¬ 
tively short integration times. In particular, as shown 
in Fig. [1] (top panel), there are top-bottom asymme¬ 
tries in the magnetic field strength that stochastically 
appear and disappear. This stochastic behavior has also 
been seen in several of the other non-dynamo cases. The 
origin of this behavior remains unclear, but it is likely 
that the asymmetry observed in Fig.[2]results from time¬ 
averaging over a small number of variations in the mag¬ 
netic held strength; averaging over many such variations 
would likely remove the asymmetry or reduce it substan¬ 
tially. 

The dynamo solution becomes more robust at larger 
radii, where the vertical depth (in terms of physical dis¬ 
tance) to which the FUV photons penetrate becomes 
larger. As a greater fraction of the vertical extent of 
the disk is ionized due to FUV photons, regions of lower 
Alfven speed overlap with high values of Am, allowing 
for the standard MRI dynamo to operate. Thus, for radii 
larger than 100 AU, we expect that the disk will be in 
the dynamo state, at least for the values of /3o explored 
here. 

As in the quasi-laminar simulations of ISimon et al.l 
(|2013all . the vertical magnetic field in the non-dynamo 
simulations can launch a magnetic wind, thus removing 
angular mo mentum through a Blandford-Payne type of 


mechanism 

Blandford & PavnellI982 

iLesur et al.ll20I3l: 

IBai & Stone 

l2013allbl; Fromane et al.l 1201311. The stress 


corresponding to this angular momentum removal lies in 
the zy component of the stress tensor and in dimension¬ 
less form is given by 


W,y 


{pVzSvy - B^By) 

Poci 


•2top 


•^bot 


(24) 


where po is the initial mid-plane gas density and Ztop and 
Zbot are the top and bottom of the magnetic wind that is 
produced as a result of the net vertical field threading the 
box. The actual values for ztop an d Zhnt, are somewha t 
arbitrary in a shearing box (e.g., ISimon et al.l l20I3all . 
and here, we take them to be the top and bottom of the 
shearing box; ztop = 4i7 and Zbot = —4iJ. This stress 
component is present in the dynamo simulations as well, 
since these simulations also contain a net vertical field. 

In the shearing-box approximation, the radially inner 
and outer sides of the box are symmetric (i.e., curva¬ 
ture is ignored). This leads to two possible types of out¬ 
flow geometry, depending on whether the outflow from 
the top and bottom sides of the box flow toward the 
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_ _ 50 AU. /Jo “ 10*. Sfui, - 0.1 g/cm* 



100 AU, fig = 10*. Sypy = 0,01 g/cm^ 
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Fig. 1.— Space-time diagram of the horizontally averaged toroidal field for two non-dynamo simulations (top and middle) and one 
dynamo simulation (bottom). For the non-dynamo runs, the toroidal field is largely constant in time at late times and either flips sign 
near the mid-plane (middle panel) or is roughly constant at all z (top panel). The dynamo simulation shows MRI dynamo oscillations for 
\z\ > 2H and a less active mid-plane region. These behaviors are representative of our simulations. 


same (“even-z” symm etry) or opposite ( “q dd-z” symme¬ 
try) radial directions (iBai &: Stonell201^ . The desired 
(physical) out flow geometry is the “eye n-z” symm etry. 
Simulations by lBai &: Stone! (|2013hD and iBail (|2013[) that 
focused on the inner region of protoplanetary disks (< 10 
AU) showed that the MRI is completely suppressed due 
to Ohmic resistiyity and ambipolar diffusion, and the 
eyen-z symmetry can be achieyed with horizontal mag¬ 
netic field flipped through a thin layer offset from the 
mid-plane. On the other hand, when the MRI is oper¬ 
ating. a s in the ideal MHD simulations of IBai fc Stonel 
(|2013aD . either the odd-z symmetry preyails (when net 
Bz is sufhciently strong), or the MRI dynamo makes the 
radial direction of the outflow change sign in time. We 
observe this latter feature in our dynamo simulations as 
well, which argues against the outflow resulting from 
a physical disk wind. Ultimately, this symmetry issue 
needs to be resolved with global simulations. Here, we 
are primarily interested in the turbulent velocity distri¬ 
bution in the disk at heights below those where a wind 
may be launched, and we need to calculate the zy stress 
only for the purposes of estimating the accretion rate 
that our disk model would support at each radius. To 
calculate the zy stress, we assume that the outflow always 
follows the physical “even-z symmetry” wind discussed 


above, which gives 




zy 


= -W 


zy 


•^top 


(25) 


Thus, the quantity of interest for our simulations will be 
the average of 2\Wzy\zt.^ and 2\WzyU,, , 


IfU.ylavg = i + 2|IU.,|.,„J , (26) 

where the factors of 2 result from the assumption made 
in equation (j25|) . and we take the average since it is not 
guaranteed that the absolute values of the zy stress are 
the same at the top and bottom of the domain. 

With this information in hand, we estimate the ac¬ 
cretion rate by assuming a steady state disk and equat¬ 
ing the stress to the accretion rate throu gh the angu¬ 
lar mo mentum evolution equation (see, e.g., Simon et al.l 
I2013al) . 


M = 


n 


(27) 


where here, R is the radial distance of the shearing box. 
This equation is only approximate since it was derived 
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Fig. 2. — Time and horizontally averaged vertical profile of the xy component of the Maxwell stress for the three representative simulations 
of Fig. ^ The solid lines are {—BxBy), where the brackets denote a horizontal average and the bar denotes a time average. The dashed 
lines are the Maxwell stress due to turbulent fluctuations, which is calculated by subtracting the correlation between large scale fields, 
— {Bx){By), from {—BxBy). The left panel shows a non-dynamo simulation in which the Maxwell stress is essentially entirely large scale. 
The middle panel shows a non-dynamo solution in which the Maxwell stress is small scale near the mid-plane (in this case due to magnetic 
reconnection of a current sheet) but is otherwise large scale. The right panel shows a dynamo simulation in which there is significant small 
scale stress due to turbulence at both large \z\ (though, large scale stress still dominates here) and in the mid-plane region. 


by calculating the mass accretion rate due to the xy 
stress alone and the zy stress alone and then adding 
the two rates. Furthermore, since the calculation of 
the wind stress should really be done in a global sim¬ 
ulation as discussed above, the contribution of \Wzy\^^^ 

to M should be taken with considerable caution. With 
these uncertainties in mind, we use Equation (|27|1 to 
calculate an order of magnitude estimate for M, which 
we list in Table [5] The accretion rates fall within the 
range lO“®-lO“^M 0 /yr, in general agreement with ob- 
servational constraints for typical T Tauri stars (e.g., 
iGullbring et al.l 119981 : iHartman n et al.l 11998). For the 
specific case of HD 16,8296. iMendigutia et al.l (|201.8D cal¬ 
culated a mass accretion rate M r; 4.5 x lO“^M 0 /yr, 
while also noting that previous results indicate a lower 
accretion rate [M ^ lO“®M 0 /yr) suggesting that the 
accretion rate has increased by an order of magnitude 
over the past 15 years. Our simulations reproduce these 
observed rates, and the runs with /3o = 10 “^ in particu¬ 
lar agree with the currently observed accretion rate for 
HD 163296. 

3.2.2. Turbulent Velocity Distribution 

We calculate the turbulent velocity distribution from 
our simulations in order to extract a turbulent broaden¬ 
ing parameter for input into LIME. Because the velocity 
distribution will cha nge with height above the mid-plane 
(|Simon et al.l 1^1 laD we calculate this distribution as a 
function of z. The total, thermal and turbulent, broad¬ 
ening can be written as 

Av{z) = Y^2/cboitzT(i?, z)lfimH + ^'turb(■2^) (28) 

where Uturb( 2 :) is the characteristic turbulent velocity as 
a function of height in a given simulation at radius R. As 
we show below, it is generally a reasonable approxima¬ 
tion to model the small-scale turbulent velocity field with 


a Maxwell-Boltzmann distribution, maintaining consis¬ 
tency between the thermal and turbulent components to 
the line width. Our basic approach is thus to fit Maxwell- 
Boltzmann distributions to the actual velocity distribu¬ 
tions, after appropriate filtering and time-averaging. 

To extract the turbulent velocities from our simula¬ 
tions, we must remove any large-scale velocity structures. 
In the shearing box approximation with orbital advec- 
tion the dominant Keplerian motion has already been 
subtracted off. However, large scale flows can be still be 
present. One example is bulk velocities that fill the entire 
domain, such as a net radial inflow or outflow. To remove 
these, we subtract the horizontally averaged velocities at 
each height. 


■Mred = v- (v)^ 


(29) 


where u^ed is the reduced velocity field and the xy sub¬ 
script on the angled brackets denote an spatial average 
in X and y. 

Another large scale flow often p resent in shearing box 
simulations are zonal flows (e.g. . iJohan sen et al.l 120091 


iSimon et al.ll2012l iSimon fc Armitagell2014[) . axisvmmet- 

ric radial variations in the gas density, pressure and az¬ 
imuthal velocity^ These flows may well be a physical 
effect (as suggested by the results of e.g., iBai fc Stonel 
|2014[) . but they occur on a radial scale that is large and 
poorly constraine d (indeed they are also s e en in global 
simulations, e.g. iDzvurkevich et al.l l2010t iFlock et 
l2011l:rUribe et al.ll2011l) . If they are observable, it is most 
likely as a spatially resolved axisymmetric perturbation 
to the density or velocity field, rather than as a compo¬ 
nent to unresolved turbulence. Accordingly, we remove 
them from the velocity field before calculating a velocity 
distribution. Since zonal flows show up in the azimuthal 

® Zonal flows are a result of geostrophic balance, where pres- 
sure gradients are b alanced by angular momentum gradients 
(IJohansen et 3111200911 . 
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Fig. 3. — Distribution of turbulent velocities for three representa¬ 
tive simulations as labelled on the individual plots. The turbulent 
veloc ities from the simulations are described in the text via Equa¬ 
tion lED and are denoted by the black curves. The red, dashed 
curves are the best fit Maxwell-Boltzmann distribution. The dis¬ 
tributions at the mid-plane (left curves in each plot) and at z = AH 
(right curves in each plot) are shown. Velocities are larger at higher 
altitudes. A Maxwell-Boltzmann distribution provides an excellent 
fit to the simulation data in the top panel and less so in the middle 
and bottom panels. The simulation in the bottom panel has the 
worst fit out of all of our simulations. Most simulations have a fit 
closer to the top or middle panels. 

velocity and are axisymmetric, we subtract their influ- 
ence by removing the y average of Vy at each x and z, 

^y.red “ 'l'y,red ~ {Vy,red)y (30) 

Having subtracted the large scale components of the flow, 
we are still left with a small scale velocity field that 
varies in space and time. To reduce this to a simple 
scalar function describing the dependence of turbulence 
on height requires several further steps. First, we ignore 
any anisotropy in the reduced velocity field, and calcu¬ 
late the magnitude of the velocity field from which we 
will calculate a distribution. 


= \/^'i.red+(^y.red)^+<red- 

Clearly, the (reduced) turbulent velocities in a disk are 
not guaranteed to be isotropic. However, we have calcu¬ 
lated a measure of the anisotropy by dividing the time- 
averaged X and z kinetic energies by the time-averaged y 
kinetic energy and taking the square root. We find that 
the different velocities differ by at most a factor of three 
within the time range over which we average. For our 
purposes, this is sufficiently close to an isotropic velocity 
field. 

Second, we assume that |z)| is statistically in a steady 
state so that it makes sense to calculate and use the time 
average of the velocity distribution. By examining the 
kinetic energy evolution for each simulation, we found 
that |u| only deviates from a statistically steady state 
for the simulations at 10 AU, where the kinetic energy 
can fluctuate by roughly an order of magnitude or more. 
This appears to be a result of sudden restructuring of 
the magnetic field within the domain. Our time aver¬ 
aging generally starts after initial transients have died 
away and ends at the end of each simulation. Since we 
employ simple arithmetic averaging, our averaged veloc¬ 
ity distributions will pick out velocities at the times when 
the kinetic energy is the largest. While this clearly rep¬ 
resents an uncertainty in our results, we have quantified 
the effect of this uncertainty on our results by comparing 
one of our radiative transfer calculations (as described in 
Section!?]) to a case where we assume zero turbulence for 
radii at 10 AU or less. We found that in terms of the 
integrated molecular line flux, the turbulence levels at 
radii < lOAU make no difference. 

Third, we examine whether the distribution of |u|, af¬ 
ter averaging, can be described via a Maxwell-Boltzmann 
distribution. Figure. |3| shows the distribution of \v\/cs 
(labelled |u|/cs in the figure for simplicity) at the mid¬ 
plane and at z = AH for three of our simulations. Each 
distribution is calculated at a number of independent 
time-slices and then time-averaged. The black curves 
are the distributions extracted from the simulations, and 
the red, dashed curves are a best fit Maxwell-Boltzmann 
distribution. The top panel shows a simulation that is 
fit well by a Maxwell-Boltzmann distribution, whereas 
the simulations in the middle and bottom panels are not 
fit as well0 The Maxwell-Boltzmann fit to the data in 
the bottom panel is the worst of all of our simulations, 
but is the only simulation that is fit this poorly; most of 
our simulations have a fit comparable to what is shown 
in the top two panels. We believe that this bad fit is a 
result of poor statistics from large temporal variations in 
the kinetic energy, as discussed above. Again, since this 
simulation is one of the 10 AU calculations, the uncer¬ 
tainties associated with this particular run will not have 
a significant impact on our results. 

Following the above procedure, we fit a Maxwell- 
Boltzmann distribution to the extracted turbulent ve¬ 
locity distribution at each height for all of our simula¬ 
tions. The characteristic turbulent velocity z^turb is then 

® At first glance, it may appear that one could obtain a better 
fit by changing the width of the distribution. However, there is 
only one parameter to fit, Uturbi which is the peak velocity of the 
distribution. 
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Fig. 4. — Characteristic turbulent velocity, defined as the peak of a Maxwell-Boltzmann distribution fitted to the velocity distribution, 
normalized by the sound speed, versus 2 : for the simulations at 10 AU. The black curves are the data, and the red, dashed curves are 
analytic fits using the piecewise function defined in Equation 13^ . The gray regions bounded by dotted lines correspond to the region of 
high diffusivity (from Ohmic diffusion, ambipolar diffusion, or both) above which FUV photons significantly enhance the ionization fraction. 
The turbulent velocities increase with \z\ in all cases, and this height dependence is steeper at the transition between the diffusion-dominated 
and FUV-dominated regions. 


the peak of the fitted Maxwell-Boltzmann distribution at 
each height. This velocity is what will be input as the 
turbulent broadening parameter in LIME. The turbulent 
velocity extracted from this fitting procedure is largely 
symmetric about the disk mid-plane, but there are small 
asymmetries. Since our LIME setup requires a perfectly 
symmetric broadening parameter, we have symmetrized 
the velocity profiles by taking the average of Uturb for 
z < 0 and z > 0 and using this average for both above 
and below the disk mid-plane. From this point on, we 
will use Uturb to denote this symmetrized velocity. 

We plot Uturb vs z in Figs. |4l|6l (black curve s). _ In 

agreement with previou s results (jSimon et al.l l2nilat 
iFromang fc Nels^ 120061 1. the turbulent velocity gener¬ 
ally increases with |z|, though in some of our simulations 
(e.g., R100-B4-FUV0.01-NG; top left panel of Fig.lHl), the 
velocity is relatively constant within the dead/damping 
zone. Within the MRI-active regions, this height depen¬ 
dence can be explained to first order by the increase in 
Alfven speed with incre asing |z|. Sinc e a oc 1//3 in MRI- 
driven turbulence fe.g.. iHawlev et anilQQSD . a ^ v\/c^, 
and since velocity fluctuations, 5v, are proportional to 
we have dv ^ ua. Simply put, larger Alfven 
speeds lead to more vigorous MRI turbulence, which en¬ 
hances turbulent velocities. 

The presence of Ohmic and/or ambipolar diffusion 
near the mid-plane damps MRI-turbulence, reducing tur¬ 
bulent velocities even more in this region. As a re¬ 
sult, the turbulent velocity profile becomes steeper at 


the transition between the FUV ionized layer and the 
diffusion-dominated mid-plane region; this is seen in all 
of our simulations and is in agreement with the results 
of iSimon et al.l (l2011all , which showed that the turbulent 
velocity at the mid-plane in simulations with an Ohmic 
dead zone is significantly smaller than the equivalent in 
a fully ionized simulation. We do not fully understand 
why the turbulent velocity profile is fiat in some simu¬ 
lations but more rounded in others, though we suspect 
that it is related to how turbulent energy is injected into 
the diffusion-dominated region from the active layers. A 
more detailed analysis of exactly why the turbulent ve¬ 
locities have the vertical profiles that they do is beyond 
the scope of this paper, but will be addressed in future 
work. 

As LIME requires an analytic function for Uturb vs z, 
we fit the following piecewise function to the curves in 
Figs. min] (red, dashed lines), 


fturb,fit(-z) = MAX 


MIN (ale^'/“^a3e^'/“") 



(32) 

where the various a coefficients are the free parameters to 
be fit. Thus, our fitting function consists of the piecewise 
combination of two Gaussians and a constant. The con¬ 
stant as is appropriate near the mid-plane for the curves 
where Uturb is relatively fiat in that region. For the other 
curves, as = 0, and this parameter is not even needed in 
the fitting procedure. As the figures show, the data is 
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Fig. 5. — Characteristic turbulent velocity, defined as the peak of a Maxwell-Boltzmann distribution fitted to the velocity distribution, 
normalized by the sound speed, versus 2 for the simulations at 30 AU. We do not include the simulation with grains included because the 
curve looks nearly identical to the case with out grains. The black curves are the data, and the red, dashed curves are analytic fits using 
the piecewise function defined in Equation 113211 . The gray regions bounded by dotted lines correspond to the region of high diffusivity 
(from Ohmic diffusion, ambipolar diffusion, or both) above which FUV photons significantly enhance the ionization fraction. The turbulent 
velocities increase with \z\ in all cases, and this height dependence is steeper at the transition between the diffusion-dominated and 
FUV-dominated regions. 


reasonably well fit by Equation (|32ll . 

4. OBSERVATIONAL PREDICTIONS 
4.1. Radiative Transfer Calculations 

In order to compare the simulations with observations, 
we use the Monte Carlo rad iative transfer code LIME 
(jBrinch fc Hogerheiidell201(ih in a post-processing fash¬ 
ion to generate simulated images of the disk projected 
onto the sky. In addition to the model parameters de¬ 
scribing the disk structure (Table 1), the simulated image 
assumes a distance of 122 pc and inclination to the line 
of sight of 44°. The assumed density and temperature 
structures are described in Section 2.1, along with the 
prescriptions to lower the CO abundance in the case of 
freeze-out in the disk interior and photodissociation on 
the disk surface. The velocity field is modeled as a cylin¬ 
drical Keplerian field. 

The turbulent linewidth is assumed to follow the func¬ 
tional form of Equation (1.121) , which reproduces the sym¬ 
metrized velocity derived from the simulations as a func¬ 
tion of height above the midplane to within a factor of 
^ 2. We use this parametric description of the turbu¬ 
lent linewidth primarily for the sake of computational 
efficiency^ In addition, the parametric form allows us 
to smoothly interpolate between the radii at which the 

LIME is capable of reading an arbitrary tabulated velocity 
field but must then perform repeated three-dimensional interpola¬ 
tions onto the Delaunay grid. 


shearing boxes sample (see description below). There are 
three significant assumptions and extrapolations we must 
make to graft the results of the shearing box simulations 
onto a continuous, two-dimensional velocity distribution 
in the disk. 

First, because the shearing boxes only sample three 
discrete radii in the disk, we must interpolate the veloc¬ 
ity as a function of radius. In order to do this, we perform 
a simple linear interpolation on each of the parameters 
oi through as from Equation (l32l) . By interpolating the 
parameters as a function of radius, we effectively inter¬ 
polate between the vertical velocity profiles displayed in 
Figs, min so that the turbulent linewidth is defined at 
all radii and heights above the midplane throughout the 
disk. Interior to the 10 AU radius of the innermost shear¬ 
ing box, we make the simplifying assumption that the 
vertical profile of turbulence is constant as a function of 
radius. This is unlikely to be true in reality, since radii 
between 1 and 10 AU are lik ely to be subje c t to an MRI- 
inactiye “dead zo ne” ('e.g.. iGammiellTQ^ iBai fc Stond 
l2013bt iBail I201§1 . but for the purposes of this initial 
investigation, in which we simulate near-future ALMA 
observations, these radii are likely to be too small to be 
spatially resolved and the assumption will have only min¬ 
imal if any effect on the high-velocity wings of the line 
profile. 

Second, there is ample evidence for vertical temper¬ 
ature gradients in protoplanetary disks, including spa¬ 
tially resolved ALMA observations of the disk around 
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Fig. 6. — Characteristic turbulent velocity, defined as the peak of a Maxwell-Boltzmann distribution fitted to the velocity distribution, 
normalized by the sound speed, versus 2 for the simulations at 100 AU. We do not include the simulation with grains included because the 
curve looks nearly identical to the case with out grains. The black curves are the data, and the red, dashed curves are analytic fits using 
the piecewise function defined in Equation The gray regions bounded by dotted lines correspond to the region of high diffusivity 

(from Ohmic diffusion, ambipolar diffusion, or both) above which FUV photons significantly enhance the ionization fraction. The turbulent 
velocities increase with \z\ in all cases, and this height dependence is steeper at the transition between the diffusion-dominated and 
FUV-dominated regions. 


HD 1 63296. on which we base our model (jBosenfeld et alJ 
l2013f) . Our shearing boxes, on the other hand, are ver¬ 
tically isothermal. There is no unique way to reconcile 
these differing structures. Here, we assume that the func¬ 
tional form of vtnTh{z)/cs, derived from isothermal runs, 
holds also when is itself a function of height. We there¬ 
fore calculate the local sound speed of the LIME model as 
a function of both radius and height above the midplane 
using the vertically varying temperature structure de¬ 
scribed in Section 2.1. At a given radius, we then use the 
symmetrized function describing the turbulent linewidth 
in units of the sound speed as a function of height above 
the midplane to scale the turbulent linewidth to the local 
sound speed. 

Finally, we need to consider what happens when the 
surface density of the disk drops below Efuv- For the 
turbulent simulations, the scaling factor 02 effectively de¬ 
scribes the height above the mid-plane at which the verti¬ 
cal turbulent velocity profile transitions from the higher- 
velocity Gaussian of the more strongly ionized upper disk 
layers to the lower-velocity Gaussian of the weakly ion¬ 
ized disk interior (see vertical lines in Fig. |6]and the ar¬ 
guments in Section [3.2.1l as to why radii larger than 100 
AU will be in the dynamo state). The vertical extent 
of the lower-velocity Gaussian should gradually taper to 
zero with increasing radius as the disk surface density 
decreases to the value of EfuVi £it which point the full 
vertical column of the disk will be maximally ionized and 
better described by the high-velocity Gaussian. To ap¬ 


proximate this effect, we gradually taper the value of 02 
beyond 100 AU so that it reaches a very small value (0.01, 
since a value of zero would cause the function to go to 
infinity) at the radius at which the surface density equals 
Sfuv (459 AU for Efuv = 0.01 gcm“^ and 308 AU for 
Efuv = O.lgcm”^). Beyond this radius, the vertical 
profile of turbulence (parameterized as a fraction of the 
local sound speed) is assumed to be constant with radius. 

4.2. Signatures of Turbulence 

In this section we perform synthetic observations of 
the turbulent disk model and examine the signatures of 
turbulence in the disk. In Section 4.2.1 we examine the 
spatial and spectral signatures of turbulence in the disk, 
both qualitatively and quantitatively, within the context 
of our fiducial model. In Section 4.2.2, we discuss the 
utility of observing multiple molecular tracers to charac¬ 
terize the vertical profile of turbulence in the disk, and 
in Section l4.2.3l we examine the effects of angular resolu¬ 
tion on the strength of the signatures of turbulence for 
the fiducial model. Throughout, we base our model on 
the properties of the bright and well-characterized disk 
around HD 163296, assuming an inclination to the line 
of sight of 44°. The observability of turbulence for any 
given set of observational and source characteristics will 
naturally depend on a host of parameters specific to the 
situation, including the specific temperature and den¬ 
sity structure of the disk being observed, the distance to 
and viewing geometry of the system (particularly inch- 
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nation), the angular resolution, spectral resolution, and 
sensitivity of the observations, and the properties of tur¬ 
bulence in the disk. Here we illustrate qualitative and 
quantitative results from a single representative choice 
of disk parameters, for the five turbulence simulations 
presented in the paper, considered against typical obser¬ 
vational parameters for a few-hour ALMA data set. 



Fig. 7. — Synthetic observations of the CO(3-2) spectra corre¬ 
sponding to a standard density and temperature model of the disk 
around HD 163296, with the turbulent velocity structure predi cted 
by the various turbulent model conditions described in Section f3.1l 
Greater turbulent velocities result in a brighter and broader line 
profile. The different model parameters are described in the leg¬ 
end. In the upper left, the scale bar indicates the magnitude of 
the 20% systematic flux uncertainty that is currently the best con¬ 
servative estimate for the accuracy of the absolute flux calibration 
with ALMA, which is comparable to the predicted difference in flux 
between the various model cases. The relative flux uncertainty be¬ 
tween channels is far smaller, and depends on the signal-to-noise 
ratio achieved by the combination of sensitivity and angular reso¬ 
lution in any particular data set. 


4.2.1. Observational Signatures of Turbulence in CO (3-2) 

We begin by generating synthetic spectra for each of 
the diffe rent /3o and Spuv conditions described in Sec¬ 
tion 13.11 Fig. [7] shows the spectral output of the LIME 
model of CO J=3-2 emission from the HD 163296 disk, 
using the fiducial temperature and density structure and 
imaged at full spectral (44 m/s) and spatial (10 AU) res¬ 
olution, for each of the five different turbulent models. 
No noise is added to the spectra, and each channel is 
summed spatially across the full 12” x 12” simulated sky 
image. For comparison, the spectrum of an otherwise 
identical disk with no turbulent linewidth is also shown. 

In general, the models with greater turbulent velocities 
yield higher fluxes and broader line peaks. The spectral 
broadening is easily understood as the redistribution of 
flux across velocity channels due to the increased tur¬ 
bulent motion of the medium. The overall increase in 
flux is best understood as an optical depth effect: the 
extra velocity gradient caused by the turbulence effec¬ 
tively spreads some of the flux from each channel into 
the neighboring velocity channels, and in a highly opti¬ 
cally thick medium there is always more CO under the 



Fig. 8. — Same as Fig. [7] but in this case the spectra have been 
normalized to demonstrate how turbulence changes the relative 
shape of the line. The systematic flux uncertainty does not affect 
determination of the relative shape of the line, which is instead 
determined by the angular resolution and sensitivity of the data. 
Diagnostics like the peak-to-trough ratio, which compares the peak 
flux of the line to the flux at line center, are useful for determining 
the amount of turbulence present in the protoplanetary disk. 

r = 1 surface that would otherwise be obscured and 
therefore unobservable. The net effect is that each chan¬ 
nel of the line is effectively spatially broadened, and the 
increase in surface area causes an increase in the flux of 
the line. Essentially, CO that would otherwise be invis¬ 
ible in the optically thick disk interior is made visible 
due to the extra velocity shear caused by the disk tur¬ 
bulence. This increase in flux is substantial; the model 
with the largest turbulent velocities exhibits a peak flux 
^22% larger than that of the model with no turbulence. 

Such a dramatic increase in flux may appear easy to 
detect, but the systematic flux uncertainty of millimeter 
data is in fact comparable to the expected flux increase, 
as is illustrated by the scale bar in the figure. This scale 
factor ultimately arises from uncertainties in the millime¬ 
ter flux models of solar system objects, which are typi¬ 
cally used to calibrate the absolute flux scale of the data 
(unlike quasars, which exhibit dramatic flux variations 
with time, solar system objects are thermal emitters with 
relatively predictable fluxes). Currently, the millimeter 
flux models of solar system objects are assumed to be 
accurate to within approximately 20%, although a goal 
is to improve these models to achieve lower systematic 
uncertainties in future generations of ALMA data. The 
systematic flux uncertainty is an uncertainty on the ab¬ 
solute flux scale of the data, and manifests itself as a 
scale factor by which the entire spectrum is multiplied. 
It therefore afflicts only the absolute value of the flux, 
rather than the relative uncertainty between different 
channels in the spectrum. The channel-to-channel rms 
noise level is instead set by the sensitivity and angular 
resolution of the data, which will differ from one data 
set to the next depending on the integration time, chan¬ 
nel spacing, and the baseline lengths in the array. As 
these factors are more readily controlled by the observer, 
measurement of the line shape rather than its overall 
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brightness is a more promising avenue for revealing the 
presence of turbulence in the disk. 


Peak-to-Trough Ratio ^1.75 

No turbulence 


/?o=l0S g/cm^ 


/?o=10" , g/cm^ 


/3o = 10 ^ , g/cm^ 


/3o =10^ . SjT^= 0 ,l g/cm^ 

T-20% T-10% Nominal T T+10% T+20% 


Fig. 9. — The peak-to-trough ratio, i.e., the ratio of the peak 
integrated line flux to the integrated flux at line center, for a range 
of models and assumed temperatures. The five turbulent model 
conditions, each of which exhibits a different turbulent linewidth, 
are arranged along the y-axis in order of most to least turbulent 
in the CO(3-2) line. The relative amount of turbulence is inferred 
from the spectral brightening and broadening evident in Fig. [3 
and the ordering of models is the same as in that figure. Along 
the x-axis, we vary the temperature structure of the disk relative 
to the nominal model conditions (in the center of the plot), up to 
a maximum of ±20% of the nominal temperature structure. The 
range of temperatures is chosen to be consistent with the uncer¬ 
tainty on the temperature structure due to the typical systematic 
flux uncertainty of 20% at submillimeter wavelengths. The peak- 
to-trough ratio is substantially more sensitive to turbulence than 
to temperature over the range of temperatures representative of 
the systematic uncertainty. 



Fig. |5] shows the same model line profiles as in Fig. El 
but normalized by their peak brightness to emphasize 
the line shape. It is evident from this figure that the line 
shape varies substantially with the amount of turbulence 
in the model, primarily in the degree to which the line 
center is “filled in” by redistributed flux from the peaks: 
the more turbulence, the more flux is redistributed be¬ 
tween channels, and the less flux contrast there is be¬ 
tween peaks and center. The “peak-to-trough” ratio is 
therefore a reasonable diagnostic of the amount of turbu¬ 
lence in the disk, independent of the absolute flux scale. 
The peak-to-trough ratio for the most turbulent model 
is 1.47, compared to a value of 1.69 for the model with 
no turbulence. This difference is detectable as long as 
the relative flux uncertainty, or channel-to-channel noise 
in the data, is limited to a few percent. Even in early 
science operations, ALMA observations routinely achieve 
this flux uncertainty on bright CO lines in nearby pro¬ 
toplanetary disks in as little as a few hours or minutes, 
depending on the details of the configuration and spec¬ 
tral resolution of the data. 

In Fig. El we examine the peak-to-trough ratio as a 
tracer of turbulence, even in the case of significant un¬ 
certainty in the disk temperature. The systematic flux 
uncertainty of 20% induces a temperature uncertainty 
that is roughly proportional, since in the Rayleigh-Jeans 
regime the line flux should be approximately linearly pro¬ 


portional to the temperature of the optically thick emit¬ 
ting surface. We therefore calculate the peak-to-trough 
ratio (the peak flux of the line divided by the flux at line 
center) for all of the different models. The models are 
arranged along the y-axis of the plot from most turbu¬ 
lent to least turbulent (according to the ordering derived 
from Figures El and E]), while varying the temperature of 
the nominal disk model by up to ±20% to simulate the 
effects of the systematic uncertainty. The range of peak- 
to-trough ratios for a given model is small compared to 
the difference in peak-to-trough ratio for a given assumed 
underlying temperature distribution and different turbu¬ 
lent conditions. While there is some small amount of 
degeneracy, the effects of turbulence dominate over the 
effects of the assumed temperature structure over the 
range of temperatures consistent with the typical sys¬ 
tematic flux uncertainty of submillimeter data. 

The effect of increasing the temperature of the disk 
is spectrally similar to the effect of increasing turbulence 
(both tend to brighten and broaden the line), however an 
increase in temperature leads primarily to an increase in 
surface brightness of the optically thick disk while an in¬ 
crease in turbulent linewidth manifests itself primarily as 
an increase in surface area. This is illustrated in Fig. 1101 
which demonstrates how the morphology of the channel 
maps changes as a function of temperature and turbu¬ 
lent linewidth. The figure displays a selection of channel 
maps from different regions of the spectrum - the line 
center (left), the line peak (center), and the line wings 
(right) - for the most turbulent model (bottom row), a 
model with no turbulence but 20% higher temperature 
than the fiducial model (middle row), and a model with 
no turbulence and fiducial temperature (top row). Line 
center refers to the channel corresponding to the systemic 
velocity, i.e., the channel with no redshift or blueshift 
along the line of sight relative to the motion of the star. 
Line peak refers to velocities near that at which the maxi¬ 
mum integrated flux occurs (in the case of our simulated 
images of HD 163296, these channels occur at roughly 
±1.3km/s relative to the systemic velocity). The line 
wings refer to the channels at velocity extremes, showing 
the most redshifted and blueshifted parts of the disk rela¬ 
tive to the systemic velocity. For Fig. 10, we have chosen 
to illustrate channels that are at roughly ±2.6 km/s rel¬ 
ative to the systemic velocity. Each panel of Fig. 13 
plots a single 44 m/s wide channel. Both at line center 
and in the line wings, the morphology of the emission in 
each channel is spatially broadened by the presence of 
turbulence in the disk. This can be understood as part 
of the process of redistribution of flux between channels. 
The morphology changes with velocity across the line, 
moving spatially from the peanut-shaped central chan¬ 
nels that trace the disk major axis to the ring-shaped 
channels that bracket the major axis near the line peak. 
When turbulence spreads the flux in the spectral domain, 
it allows emission to “bleed” between nearby channels, 
effectively broadening the range of spatial locations that 
are represented within each channel. The result is the 
relatively broader morphology illustrated by Fig. 1101 in 
which the models illustrated in the top and bottom rows 
have an identical temperature structure but differ in the 
amount of turbulence that they exhibit. 

Interferometric molecular line data is fundamentally 
three-dimensional in nature, and any fit to the data 
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Fig. 10.— A selection of channel maps generated for simulated images of the turbulent disk models at high (0.08 arcsec) angular 
resolution. The left column illustrates the morphology of CO(3-2) emission at line center, i.e., the systematic velocity. The center column 
illustrates the morphology near the line peak, and the right column the line wings. The top row shows the fiducial disk model with no 
turbulence, the middle row shows the same model with no turbulence but 20% higher temperature (the maximum allowed by the systematic 
flux uncertainty), and the bottom row shows the model with the strongest turbulence in the upper disk layers probed by the CO(3-2) line 
(/3o = 10'*, Spuv = 0-1 g cm~^). The primary signature of increased turbulence in channel maps is the spatial broadening of emission, 
whereas increased temperature only makes the emission brighter. 


would naturally take into account the full position- 
position-velocity cube, including constraints both from 
the morphology of the line at each velocity and the over¬ 
all shape of the spectrum. The simulated images of the 
model illustrate that both the spectral line shape (in¬ 
dependent of the uncertain absolute flux scale) and the 
morphology of the emission within each channel provide 
constraints on the amount of turbulence in the disk, and 
that the turbulent broadening can be distinguished from 
temperature broadening given sufficient relative flux un¬ 
certainty and spatial resolution. 

4.2.2. The Role of Different Molecular Traeers 


Section 4.2.1 illustrated signatures exclusively from the 
CO(3-2) line, which is one of the most commonly ob¬ 
served tracers of molecular gas in circumstellar disks. 
The advantages of observing CO(3-2) are that it is 
bright and occurs in a relatively optically thin window 
of Earth’s atmosphere; it is therefore a powerful and 
commonly-observed tracer of the gas component of cir¬ 
cumstellar disks. However, in the context of understand¬ 
ing the turbulent structure of protoplanetary disks, its 
use is limited. It is extremely optically thick and probes 
primarily the upper surface layers of the disk, with the 
emission originating from a height of 3-5 scale heights 
above the midplane in our fiducial models. A judicious 
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Fig. 11.— Cartoon showing the layered structure of the disk and illustrating how the optical depth determines the height above the 
midplane from which different molecular tracers originate. In the uppermost surface of the disk, CO is photodissociated by energetic stellar 
radiation and cosmic rays. In the cold midplane, CO freezes out onto solids and therefore does not emit the gas-phase lines. CO exists 
only in the warm molecular layer. Very optically thick transitions in the CO rotational ladder emit from a region very close to the surface 
of the warm molecular layer. CO isotopologues, which have a lower optical depth, have their r = 1 emitting surface deeper in the disk, 
where the disk is colder and less turbulent. In our model of the HD 163296 disk, the densities are high enough that ^®CO is optically thick, 
while C^®0 remains optically thin throughout the vertical column of the disk. 


choice of molecular line tracers can provide insight into 
the three-dimensional structure of the disk, allowing the 
observer to probe much closer to the planet-forming re¬ 
gions near the disk midplane. 

Carbon monoxide is the molecule of choice for most 
disk observations due to its status as the second most 
abundant molecule in protoplanetary disks (the most 
abundant molecule, H 2 , is virtually invisible due to its 
lack of a dipole moment). Observations of CO isotopo¬ 
logues, which have lower optical depths than ^^C^®0, 
allow us to probe deeper into the warm molecular layer, 
as illustrated by the cartoon in Fig. [TT] Protoplanetary 
disks exhibit a vertically and radially stratified struc¬ 
ture, including a hot, atomic upper skin in which CO 
and other molecules are photodissociated; a cold, densely 
shielded midplane within which volatile species freeze 
out of the gas phase onto the surfaces of dust grains; 
and a warm molecular layer in which the conditions are 
neither to o energetic nor too cold for the formation of 
CO (e.g. lAikawa et al.l 119961 : lAikawa fc Nomural 120061 : 
iGorti fc Hollenbachll2008[l . Efforts at modeling the den¬ 
sity structure of disks reveal that in a typical bright pro¬ 
toplanetary disk, both ^^C^®0 (the most common iso¬ 
topoiogue, generally abbreviated as just “CO”) and its 
next-most-abundant isotopoiogue ^^C^®0 (generally ab¬ 
breviated “^^CO), are optically thick, although ^^CO 
is more optically thin and therefore origin ates from a 
surface that is deeper within the disk fe.g.. iPanic et al.l 
[ 2 OOI). 12^180 (generally abbreviated C^®0) is the 

next-most-abundant isotopoiogue, and it is typically the 
brightest optically thin tracer, dominated by the dense 
regions close to the midplane but at temperatures too 
high for CO freeze-out to have occurred. It is therefore 


possible to build a vertical profile of turbulence through¬ 
out the disk by examining how the peak-to-trough ratio 
varies for different molecular tracers. 

This approach is illustrated by Fig. [121 The top row 
shows simulated spectra of three isotopologues in terms 
of the absolute flux, while the bottom row shows normal¬ 
ized flux. The left column shows CO(3-2), which probes 
the uppermost surface of the warm molecular layer, while 
the center shows ^^CO, which probes the middle of the 
warm molecular layer, and the right column shows C^®0, 
which probes the innermost regions of the warm molec¬ 
ular layer close to the CO freeze-out zone in the mid¬ 
plane. The model with greatest turbulence (/3 = 10^ 
Spuv = 0.1 g cm“2) is represented by the solid blue 
line, the same model with no turbulence is represented 
by the dashed green line, and a model with a constant 
turbulence level of 0.46Cs, chosen to match the normal¬ 
ized CO(3-2) line profile of the most turbulent model, 
is represented by the dotted red line. It is evident that 
the peak-to-trough ratio increases from the least to most 
optically thick line tracer, reflecting the increase in tur¬ 
bulent linewidth with height above the midplane. It is 
also clear from the figure that the differences in spec¬ 
tral shape between the MHD models in which turbu¬ 
lence varies as a function of height above the midplane, 
and the constant-turbulence models in which turbulence 
is a constant fraction of the sound speed (and therefore 
depends only on temperature), is more subtle, especially 
in the context of the relativ ely large ALMA s ystem atic 
flux uncertainty. A paper bv iHorne fc Mar^ ( 19861) de¬ 
scribes the relative lack of sensitivity of optically thin line 
spectra to levels of turbulence in accretion disks, and em¬ 
phasizes the importance of the spatial domain, which can 
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Fig. 12.— Different molecular line tracers with varying optical depth can probe the vertical profile of temperature through the disk. The 
top row shows the absolute flux as a function of velocity predicted by the p = lO'^ SpuY = 0-1 g cm“^ model (solid blue line) compared 
with a model with no turbulence (dashed green line), and a model with a constant turbulent linewidth as a function of sound speed, 0.46 Cg 
throughout the disk, chosen to match the normalized CO(3-2) line profile of the most turbulent model (red dotted line). The bottom row 
shows the same spectra, normalized to emphasize the peak-to-trough ratio. The left column shows the predicted CO J=3-2 spectra (highest 
optical depth), the middle column shows the J=3-2 spectra (intermediate optical depth), and the right column shows the J=3-2 

spectra (lowest optical depth). The difference in the peak-to-trough ratio for the turbulent vs. non-turbulent models is greatest for the 
high optical depth tracer that shows the more MRI-active surface of the disk, and smallest for the low optical depth tracer that shows the 
conditions closer to the midplane. This reflects the turbulent velocity gradient that incr ease s with distance from the midplane, depicted in 
Fig. El as well as the effect of the layered disk structure, depicted by the cartoon in Fig. EU The difference between a constant-turbulence 
model and a variable-turbulence model is more subtle, especially in the context of the normalized line profiles used to circumvent the 
ALMA systematic flux uncertainty. 


provide measurable signatures of turbulence even in the 
absence of differences in the unresolved line profile. 

While we do not simulate them here, there are other 
molecular line tracers that are advantageous for probing 
the three-dimensional structure of turbulence in the disk. 
Heavier molecules, like CS, exhibit a naturally lower ther¬ 
mal broadening than CO due to their reduced molecular 
weight, and therefore turbulent broadening is easier to 
disentangle from thermal broadening (see Guilloteau et 
al. 2012). The disadvantage of such a rare species is that 
in an optically thin line temperature and density become 
significantly degenerate, which adds another degeneracy 
to the turbulent linewidth determination; furthermore, 
the chemistry of CS and its spatial distribution vertically 
and radially through the disk is less well understood. An¬ 
other promising tracer is DCO"*", which is destroyed by 
CO and therefore becomes abundant in regions where CO 


is frozen out in the cold midplane. DCO'*' is therefore a 
promising tracer for placing constraints on turbulent line 
widths in the midplane of protoplanetary disks. 

4.2.3. The Role of Resolution 

Observational characterization of turbulence in disks 
requires both high angular resolution and high spectral 
resolution. The spectral resolution should of course be 
smaller than the expected turbulent linewidth to ensure a 
detection (see, e.g., Isella et al. 2007); however, angular 
resolution is equally important to aid in disentangling 
the effects of turbulence from those of temperature. 

Figure IT5I illustrates the importance of angular resolu¬ 
tion in characterizing the turbulent linewidth. For our 
fiducial model (based on the observational parameters of 
HD 163296), we examine three slices in RA across the 
disk in the central channel (located at the systemic ve- 
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locity), at projected linear separations of 100, 200, and 
300 AU from the central star. The locations of the slices 
are illustrated in the leftmost panel, superimposed over 
a disk model with no turbulence (taken from the top left 
panel of Fig. uni). The right three panels show flux as 
a function of RA offset across each slice. In each panel, 
a model with no turbulence (solid line) is compared to 
the most turbulent model from our simulations (/3 = 10^, 
Spuv = 0.1 g cm“^). The absolute values of the flux as a 
function of position are highly dependent on the details of 
the model (viewing geometry, turbulent linewidth, etc.) 
and the parameters of the observation (particularly spec¬ 
tral resolution); however the trend of broader emission 
with greater turbulence is clearly evident. In the bot¬ 
tom row of the figure we highlight detectability with a 
realistic simulation of ALMA observations of the model, 
assuming three hours of observation at 0"2 resolution 
with the sensitivity anticipated for Cycle 3. While spa¬ 
tial filtering does somewhat affect the overall flux levels, 
the trend of increasing spatial width with turbulence is 
readily evident in realistic simulated observations. 

The different radii show an important trend that em¬ 
phasizes the necessity of high angular resolution when 
attempting to characterize turbulence at radii close to 
the central star. From the first to the third slice, the dis¬ 
tance from the star triples, and the absolute width of the 
line approximately doubles (as expected for the spatial 
separation of isovelocity contours of a disk inclined at 44° 
to the line of sight). The difference in spatial width be¬ 
tween the turbulent and non-turbulent models over the 
same range of radii increases more slowly, however, from 
^0"4 at 100 AU to ^Ol'O at 300 AU. This is also to be 
expected, since the turbulent linewdith is expected to be 
proportional to the sound speed, which is proportional 
to the square root of the temperature, and the tempera¬ 
ture decreases roughly as ^Iffr. For a constant turbulent 
linewidth as a fraction of the sound speed, therefore, the 
velocity width of the line should change quite slowly with 
distance from the star, approximately as this will 

result in a more slowly-changing physical width of the 
emission as viewed in projection on the sky. Of course, 
this scaling neglects important properties of tur¬ 

bulence that are included in our models, like the chang¬ 
ing ionization fraction and UV penetration depth as a 
function of radius, and the more shallow radial temper¬ 
ature dependence due to the flaring of the disk. How¬ 
ever, the direction of the trend holds: the turbulent 
linewidth should change relatively slowly with position 
in the disk, indicating that the physical width of the line 
should change relatively slowly with distance from the 
star. 

This latter trend is good news for studies of turbu¬ 
lent linewidth in the inner disk: as in Fig. [T31 the to¬ 
tal spatial width of the line will decrease rapidly as one 
moves towards the central star, but the spatial difference 
in linewidth between the turbulent and non-turbulent 
models will fall off more gradually. Due to the relatively 
small emitting area of the inner disk, the overall spec¬ 
tral shape of the line emission is dominated by the larger 
surface area of the outer disk. High angular resolution is 
therefore extremely important for detecting turbulence 
in the inner disk, primarily due to its effects on the spa¬ 
tial width of the line emission. The exact angular res¬ 


olution needed to resolve turbulence will depend on the 
details of the viewing geometry of the target, the ex¬ 
pected amount of turbulence in the disk, etc., but in our 
fiducial highest-turbulence model the spatial turbulent 
broadening at 100 AU requires an angular resolution of 
<0.4". 

5. SUMMARY AND CONCLUSIONS 

Spatially resolved ALMA observations of protoplan¬ 
etary disks have the potential to yield constraints on 
turbulence that are more direct than those available in 
other disk-accreting systems. Maximizing that poten¬ 
tial, however - and determining empirically the nature 
of protoplanetary disk turbulence - requires first under¬ 
standing in detail the observational characteristics of dif¬ 
ferent models of disk turbulence. As a first step toward 
that goal, we have presented the results of local simula¬ 
tions of magnetohydrodynamic turbulence in the outer 
disk, incorporating the effects of ambipolar and Ohmic 
diffusion. The simulations were tailored to match a rep¬ 
resentative model for the disk around HD 163296, and 
considered radii between 10 AU and 100 AU. We de¬ 
veloped and validated a methodology for extracting the 
turbulent velocity distribution from each simulation, and 
incorporating it into a radiative transfer calculation that 
was used to generate synthetic observables. 

Our work should be considered as preliminary in sev¬ 
eral regards. Most obviously, although the numerical 
simulations include what we believe to be the most essen¬ 
tial physics operating on these scales, they have known 
limitations. The simulations are local, isothermal, ne¬ 
glect the Hall effedE3, and are limited in terms of the 
range of parameters studied and length of integration. 
None of these limitations represent a fundamental obsta¬ 
cle, though relaxing them all (in the absence of algorith¬ 
mic improvements) would require computing resources 
well beyond what is currently available. We consider the 
simplified thermodynamics to be the most serious limi¬ 
tation. The assumption of an isothermal disk is not only 
inconsistent with the radiative transfer calculations but 
may also affect the propagation of waves through the at¬ 
mosphere, as wave propagatio n is known to de pend on 
the vertical thermal structure (jBate et al.ll2002h . 

In this analysis we have identified several promising 
tracers of turbulence and explored in a qualitative way 
the major known degeneracy between temperature and 
turbulence. The observers ability to constrain turbulence 
given degeneracies with other parameters depends upon 
the details of the disk and the observational parameters 
of a given data set, for example, the size, temperature, 
and viewing geometry of the disk, and the angular reso¬ 
lution and signal-to-noise ratio of the observations. The 
question of the observers ability to place strong statisti¬ 
cal constraints on the amount of turbulence in the disk 
given these degeneracies is an important one, but not 
straightforward to answer in a general way. A forthcom¬ 
ing publication (Flaherty et al. in prep) will present a 
thorough statistical analysis of turbulence measurements 
using ALMA observations of the disk around HD 163296, 
including an exploration of the degeneracies with related 

Preliminary calculations of the turbu lent velocity structure 
that include the Hall effect, bv IBail (I2014bll . suggest that this may 
not substantially impact the results in the outer disk. 
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Fig. 13.— A sketch of the role of angul ar r esolution in the detection of turbulence. The left panel shows the central channel of a model 
with no turbulence (top left panel of Fig. [TO] above), indicating cuts across the surface brightness of the disk at projected distances of 100, 
200, and 300 AU from the central star. The three panels to the right show the flux per pixel in the image as a function of offset along 
the RA cuts indicated in the leftmost panel, at distances of 100, 200, and 300 AU, respectively. The non-turbulent model (solid line) is 
spatially narrower than the turbulent model (/3o = 10 '*,Sfuv = 0.1 g cm~^; dotted line). The top row is for the non-spatially-filtered 
full-resolution simulated image generated with the Lime radiative transfer code, while the bottom row shows the same models (including 
appropriate noise) after a three-hour Cycle 3 ALMA observation in the C36-4 configuration at 0”2 angular resolution and 44 m/s velocity 
resolution. The numerical difference in width of the turbulent model relative to the non-turbulent model depends on the chosen spectral 
channel and the details of the parameters chosen (particularly inclination, spectral resolution, and turbulent linewidth); this figure serves 
merely to illustrate the importance of spatial resolution in distinguishing turbulent from non-turbulent disk models. 


parameters and robust statistical characterization of the 
degree to which turbulence can be measured given these 
degeneracies. 

Our primary conclusions are as follows: 

1. Magnetohydrodynamic turbulence in the observ¬ 
able outer regions of the HD 163296 disk is sub¬ 
stantially modified by ambipolar diffusion. The 
MRI can yield accretion rates of the same order 
of magnitude as those observed if the disk surface 
is ionized by FUV photons, and a net vertical mag¬ 
netic field is present. Weaker MRI turbulence and 
lower accretion rates are possible, but would be in¬ 
consistent (within an MRI-dominated disk model) 
with the idea that the inner disk is sustained from 
large scale accretion. 

2. The characteristic turbulent velocity generally in¬ 
creases with the vertical height above the disk 
mid-plane, \z\. The slope of this velocity gra¬ 
dient changes at the transition to fully active 
MRI turbulence in the FUV ionization layer. In 
units of the sound speed, this velocity ranges from 
■Cturb 0.01 — O.lCs in the disk mid-plane to 
f^turb ^ 0.1 — 0.4cs at the bottom edge of the ac¬ 
tive region. Within the active region, the turbu¬ 
lent velocity can reach the sound speed, Uturb Cg. 
This height dependence for the turbulent velocity 
appears to be a generic property o f MRI-driven 
turbulence (iFrornang fc Nelsonll2006t iSimon et al.l 

[2FviTa [^[2 oTh^ 

3. The outcome of the MRI in the ambipolar dom¬ 


inated outer disk yields two classes of solution, 
a “dynamo” solution in which the magnetic field 
structure in the FUV-ionized layer displays peri¬ 
odic reversals (as in the standard picture of the 
MRI in ideal MHD, e.g., ISimon et al.ll201^ . and 
a “non-dynamo” solution in which the magnetic 
field structure is relatively constant in time and 
the Maxwell stress has a substantial large scale 
component. The dynamo solution appears to be 
preferred at larger disk radii. The magnitude and 
vertical gradient of turbulent velocity, however, is 
essentially independent of the solution type. 

4. The simplest observational diagnostic of the level 
of turbulence is the spatially integrated ratio of 
the peak flux to the line center flux (the “peak- 
to-trough ratio”). For the CO(3-2) line, we find 
readily detectable variations in this ratio between 
our non-turbulent and most-turbulent models. In 
contrast to the absolute line flux (which is also a 
function of the turbulent velocity structure in the 
disk), the dependence of the peak-to-trough ratio 
on turbulence is only mildly degenerate with tem¬ 
perature. Moderate temperature uncertainties, re¬ 
sulting from systematic flux uncertainties at mil¬ 
limeter wavelengths, do not prevent a measurement 
of the turbulent broadening component. 

5. Spatially resolved disk observations yield addi¬ 
tional diagnostics of disk turbulence. Most obvi¬ 
ously, it may be possible to directly observe large 
scale structures that develop spontaneously within 
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different classes of turbulent flows (e.g. zonal flows 
in MRI turbulence, or vortices in some alternative 
models). Small scale turbulence manifests spatially 
as a broader distribution of flux in a given veloc¬ 
ity channel (assuming spectral resolution less than 
the turbulent linewidth), while higher temperature 
(within the range consistent with a 20% systematic 
flux uncertainty) primarily brightens the emission 
in a given channel. This effect is most pronounced 
at the line center and peaks, and has a smaller 
effect on the line wings. Due to the relatively shal¬ 
low dependence of turbulence on radius, the spatial 
width of the line due to Keplerian velocity shear de¬ 
creases rapidly, while the spatial width due to tur¬ 
bulence decreases more slowly. High angular reso¬ 
lution is therefore a potent probe of turbulence in 
regions close to the central star, whose relatively 
small emitting areas make spectral probes of tur¬ 
bulence relatively insensitive. 

6 . Different molecular line tracers that vary in optical 
depth trace different heights in the disk. Analy¬ 
sis of an ensemble of line tracers therefore allows 
us to probe the vertical distribution of turbulent 
linewidth as a function of height above the mid¬ 
plane. Radiative transfer calculations are needed 
to determine which lines provide the best discrim¬ 
inants between different turbulence models. 

The known properties of the HD 163296 disk, to- 
gether with existin g results from SMA observations 
(jHiighes et al.ll2niTl l. suggest that near-term ALMA ob¬ 
servations will be able to constrain protoplanetary disk 
turbulence through all three of the above diagnostics 
- the peak-to-trough ratio, spatially resolved channel 
maps, and distinct optical depth tracers. Our results 
suggest that the MRI, at least, displays distinctive ob¬ 
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